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Abstract. We analyse the time-dependence of currents in a ID Bose gas in an optical 
lattice. For a ID system, the stability of currents induced by accelerating the lattice 
exhibits a broad crossover as a function of the magnitude of the acceleration, and the 
strength of the inter-particle interactions. This differs markedly from mean-field results 
in higher dimensions. Using the infinite Time Evolving Block Decimation algorithm, we 
characterise this crossover by making quantitative predictions for the time-dependent 
behaviour of the currents and their decay rate. We also compute the time-dependence 
of quasi-condensate fractions which can be measured directly in experiments. We 
compare our results to calculations based on phase-slip methods, finding agreement 
with the scaling as the particle density increases, but with significant deviations near 
unit filling. 
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1. Introduction 

Exciting progress in experiments with cold atoms in optical lattices [Ij has not only paved 
the way for study of the quantum phases associated with strongly interacting many-body 
systems [21 13], but also for study of non-equihbrium dynamics in such systems. This 
is particularly true for transport properties, where the long coherence times associated 
with the experiments make it possible to gain new insight into phenomena such as spin- 
charge separation [H El E], and currents in the presence of impurities and junctions 
O El [3 Ho] by studying them in a new environment. Particularly in the case of ID 
systems, this connection to non-equilibrium dynamics is further strengthened by the 
recent development of numerical methods based on matrix product states [11] , including 
the time evolving block decimation (TEBD) algorithm [121 [13] and the related time- 
dependent density matrix renormalisation group (t-DMRG) methods [HI |T5]. These 
make possible the quantitative prediction of time-dependent dynamics for size scales 
that are typical in experiments, as well as the identification of parameter regimes in 
which specific phenomena can be observed. 

A key characteristic for the transport of bosons in a lattice is the stability of 
currents in the presence of interactions which can arise in regimes where the system is 
superfluid. However, this stability depends not only on the current and the interactions 
[T6l [T71 flSl [19] , but also on the dimensionality of the system. This was demonstrated 
in recent experiments by observing damping of the centre of mass motion for bosons 
oscillating in an harmonic trap [20]. For sufficiently small initial displacements, the 
motion was found to be stable up to a critical lattice depth in higher dimensions, 
whilst for bosons confined to move along one dimension, decay of the oscillatory motion 
was observed for arbitrarily small depth. This was explained by the increased role of 
quantum fluctuations in the ID system [211[22l[23l[2ll[25l[26l[27]. A similar situation 
occurs in the case where a homogeneous current is created, e.g., by starting in the lowest 
energy state, and then accelerating the lattice [281 [291 EOl ED E2] (as shown in figure [1]). 
A mean-field stability diagram as a function of the magnitude of the acceleration and 
the strength of inter-particle interactions was computed by Altman et al. [331 El], 
predicting a sharp transition between stable and unstable regions. For bosons confined 
in 3D, such transitions were observed recently [31]. However, the same experiments 
indicated strong deviations from this behaviour with a ID gas, for which a crossover 
was observed between the two regimes. 




Figure 1. Bosons in the ground state configuration in an optical lattice are 
instantaneously accelerated to a mean quasi-momentum in the lowest Bloch band. 
The stability of the resulting current will depend on the value of this quasi-momentum 
and the inter-particle interactions. 
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Here we perform a detailed numerical analysis of the decay of currents for a Bose 
gas moving uniformly in ID, making quantitative predictions for the time-dependence of 
currents that can be measured in an experiment. This analysis is made possible via the 
use of the infinite TEBD (iTEBD) algorithm |35], which until now has primarily been 
applied to study the ground states of homogeneous systems. Here it is used to investigate 
the time-dependence of homogeneous current flow, allowing predictions to be made over 
longer timescales without boundary effects. We begin from the ground state of the Bose- 
Hubbard model [2], which describes the lattice gas, and consider an acceleration applied 
in order to produce a finite current. As a function of the interaction strength in the gas 
and the magnitude of the acceleration, we observe a broad crossover between regions 
of stable and unstable current, which we characterise via the variation in the rate of 
decay of the current for different parameters. We find that the crossover region occurs 
at significantly smaller values of interaction strength and initial acceleration than are 
predicted by mean-field methods. We also compare our decay rates to those obtained for 
systems with large filling factors via phase-slip calculations by Polkovnikov et al. |34] . 

The rest of this paper is organised as follows: In section |2] we discuss the Bose- 
Hubbard model for the system we study, and the methods we use to characterise the 
currents in the system. In section [3] we describe the use of the iTEBD algorithm and 
the procedure we use in our calculations. In section H] we present a quantitative analysis 
of the decay of currents within the system, as well as the dependence of the decay on 
the initial acceleration and the inter-particle interactions. The scaling properties of 
the decay rates are investigated in section |5l and compared to calculations based on 
phase-slip methods, and in section [6] we give a summary and outlook for this work. 



Appendix A provides detailed information about translationally invariant block sizes 



within the iTEBD algorithm. 



2. Moving ID bosons in an optical lattice 

In this section we introduce the system of bosons moving in ID in an optical lattice. 
We first introduce the Hamiltonian describing the system in 12.11 before discussing the 
current and general expectations for its stability during time evolution in 12.21 In order 
to make stronger connections to quantities that can be measured in an experiment, we 
introduce the quasi-condensate fraction in 12.31 the decay of which is strongly related to 
the decay of the current. 



2.1. Bose- Hubbard model 

The dynamics of bosons in the lowest band of an optical lattice is described by the 
Bose- Hubbard model [H 121 [3] with Hamiltonian {h= 1): 

H = -jJ2brb^ + ^Y.^^{n,-l). (1) 
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Here, hi {hi) are the bosonic annihilation (creation) operators at lattice site i, hi = h\hi^ 
U is the on-site interaction energy shift, and J is the tunnelling amplitude between 
neighbouring sites, with ^^j^) denoting a sum over neighbouring lattice sites. This 
model is valid in the regime where Uh, J <^ ut, where ut is the energy separation 
between the lowest two Bloch bands, and we denote the filling factor of particles on 
M lattice sites as = N/M. 

For integer n, a phase transition is observed in the ground state of this model as a 
function oiu = U/J between superfiuid and insulating behaviour [36]. In one dimension 
with n = 1, the critical value of m is Mc ~ 3.37 [37]. The superfiuid (SF) ground state 
is characterised by quasi off-diagonal long range order in the single particle density 
matrix (SPDM) {blbj). For u > Uc, the system is in a Mott insulator (MI) phase, with 
exponentially decaying off-diagonal elements of the SPDM Q)\bj). In what follows below, 
we will assume that the system is prepared in the ground state of the Bose-Hubbard 
model for a particular choice of n and m, and then accelerated to produce a finite initial 
current. 



2.2. Currents in the Bose-Hubbard model 

By accelerating the lattice or applying a linear gradient potential for a short period 
of time, it is possible to create a current of atoms moving with respect to the optical 
lattice. This can be quantified via the current operator 




which appears in the continuity equation 

jk-jk-i = ^hk. (3) 

Note that the average of the current expectation values ^m(irn) /M is proportional to 
the mean group velocity calculated via the quasi-momentum distribution of the particles. 
In particular, we can define bosonic operators for the quasi-momentum modes, Og, which 
are related to hi as 6^ = Ylq e~*''"™ag/ a/M, where a is the lattice spacing, and q = 27rr/L 
for r G (— M/2, M/2] an integer, with L = aM the lattice length. We see that 
1 ^ 1 

^ J2^Jm) = J^Y^'^J sin (ga)ng, (4) 

m q 

where Uq = (ajciq) is the quasi-momentum distribution function at quasi-momentum q, 
for which the corresponding group velocity per lattice constant is given by 2Jsin(ga). 
In this study, the current will always be translationally invariant, and so we omit the 
site label and write (j) = {jk), for any k. 

If a finite current is induced to the system by accelerating the lattice so that the 
quasi-momentum distribution is shifted by an amount ka (see figure [1]), then we expect 
the resulting behaviour to depend on the superfiuidity of the gas. Even in the presence 
of interactions, a superfiuid fiow without dissipation can exist. However, if a current is 
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generated in a non-superfluid initial state, or if the flow becomes dynamically unstable, 
then the current will decay. Whether the current is stable or not is thus a function of 
both the initial state of the system and the initial mean quasi- momentum, ka |34j . 

As has been predicted theoretically in the weakly interacting Gross-Pitaevskii 
regime [HI [19] and demonstrated in experiments [171 EHl EH] , the current will be unstable 
if the initial superfluid ground state (with all atoms near ka = 0) is accelerated to 
the inverted part of the lowest Bloch band, i.e., when ka > 7i/2 (this corresponds 
to a classical instability). We thus expect stable currents only when ka < 7t/2. If 
we begin with the ground state of a system with integer filling n, then we also do 
not expect stable superfluid currents above the critical interaction Uc for the SF-MI 
phase transition. A stability diagram, interpolating between these two regimes was 
investigated by Altman et al. [33] using a Gutzwiller mean-field technique, and they 
found a sharp transition between stable and unstable regions. This sharp transition 
was observed in an experiment by Mun et al. [31] for bosons allowed to move in the 
lattice in three dimensions. However, the same experiment observed behaviour more 
characteristic of a crossover between stable and unstable currents in ID. 

In this sense, in a ID system it is not possible to characterise clear regions where the 
currents decay and do not decay, but rather the crossover must be characterised by how 
rapidly the current decays for different initial mean quasi-momenta and inter-particle 
interaction strengths. In section H] we will investigate this crossover quantitatively in 
terms of such decay rates using the iTEBD algorithm. 

2.3. The quasi- condensate fraction 

In section H] we show that the decay of the current results from two-particle scattering 
processes that broaden the quasi-momentum distribution. In 2D or 3D decay of the 
current is thus also directly linked to a decrease in the condensate fraction in the 
experiment, which was used as the key experimental observable in reference [31]. In 
ID, there is no condensate for an infinite system, however as the system is finite (of 
the order of 100 occupied lattice sites), a peak that would correspond to the condensate 
is observed in time of flight measurements of the momentum distribution. In order to 
make connection to this observable, we will calculate not only the time dependence of 
the current (j) below, but also a quasi- condensate fraction defined over a finite portion 
of the system. This fraction, Cr, can be defined as the largest eigenvalue of the reduced 
SPDM considering only a range of R sites in the system i.e. (blbj), with i,j G [/, / + R] 
for some arbitrary site I. We typically use R = 100 below, reflecting the typical occupied 
number of lattice sites in experiments. 

3. Numerical Calculations of the current decay 

In the remainder of this manuscript, we analyse the decay of a current, which is created 
by suddenly imparting a quasi-momentum ka at time t = to bosons loaded into the 
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lattice in the ground state configuration (which has zero mean quasi-momentum) . This 
resembles the experimental situation of reference [31], where a current is obtained by 
means of a moving optical lattice which is accelerated to a final quasi-momentum ka on 
a timescale long enough to ensure adiabaticity with the respect to inter-band transitions. 
Here we assume that the acceleration is fast compared with the tunnelling timescale 1/ J, 
resulting in a simple translation of the quasi-momentum distribution, (see section 14.11 
and figure [2]). This corresponds to an application of the operator 

k{ka) = Yl e-*('='^)'"' (5) 
I 

to the initial ground state. 

To compute both the initial ground state of the system, and the time dependence 
of the current once the initial mean quasi-momentum ka is imparted on the system, we 
make use of the iTEBD algorithm [35]. This algorithm is an extension of the TEBD 
algorithm [T2l [13], and makes possible the near-exact integration of the Schrodinger 
equation for ID lattice Hamiltonians in an infinite, translationally invariant system. 
These methods have not only been applied to coherent dynamics in ID, but also 
generalised to dissipative systems [38l [39l HQ], and extensions to a 2D state ansatz 
have been recently considered |llj . 

Until now, iTEBD has been primarily used to compute the ground state of 
translationally invariant systems, however here this algorithm becomes crucial in 
computing time-dependent dynamics 0- By computing the decay of the current in 
an effectively infinite system, we avoid difficulties arising due to either open boundary 
conditions, which can prevent propagation of moving bosons, and periodic boundary 
conditions, which restrict the length over which correlation functions can be computed. 
In our calculations we apply the iTEBD algorithm as described in reference [35], but 
with the modification that the length of a single block in the translationally invariant 
system is increased from two to 271 /ka. Although a matrix product state representation 
can reproduce phase relationships with any period, even with a single repeated matrix, 
we find it convenient to increase the block size so as to represent one full cycle of the 
phase when we apply the operator K{ka) to the state (see Appendix A for more details). 

Throughout our calculations, we performed convergence tests in the matrix product 
state bond-dimension x, as well as the time step for integration, and the maximum 
number of particles allowed per site in our description of the Hilbert space, d. For a 
more detailed analysis of how errors enter this method, see 



4. Analysis of Current decay 

In this section we characterise the current and its decay as a function of the initial mean 
quasi-momentum ka and the strength of interactions m, as well as the evolution of the 
characteristics of the state as a function of time. We begin in section l^TTl by analysing the 

I Time-dependent dynamics with the iTEBD algorithm have also been recently employed to study 
quantum quenches [41]. 
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initial state of the system at a time t = 0+, after the instantaneous acceleration of the 
lattice results in a mean initial quasi- momentum ka. In section [^^2] we then describe the 
dynamical changes in the quasi-momentum distribution during time evolution. Finally, 
we discuss the dependence of the decay on ka and u and provide a shaded plot of a 
stability crossover diagram in section 14.31 

4.1. Initial current 

In this subsection we characterise the system configuration at a time t = immediately 
following the initial instantaneous acceleration (corresponding to an imposition of the 
initial mean quasi-momentum ka). We consider both the quasi-momentum distribution 
Uq and the corresponding current (j). 

0.7 
0.6 
0.5 
^ 0.4 

0.3 
0.2 
0.1 


-0.1 0.1 0.2 0.3 0.4 0.5 

q/iT 

Figure 2. Normalised quasi-momentum distributions before {t = 0^) and after 
{t — 0+) a momentum ka — 0.47r is imparted on the particles. The solid lines show 
the distributions for an inter-particle interaction of m = 1 (SF regime), the dotted lines 
are for distributions with u = 4 (MI regime). The numerical parameters are x — 100) 
d = 6. 

Figure |2] shows Ug before and after the instantaneous application of a momentum 
shift of ka = OAn at t = for a situation where the lattice filling is n = 1. At t = 0~ 
(that is, before the shift is applied) the system is in a SF or MI configuration, with 
u = 1 and 4, respectively. At time t = 0~ and for u = 1, rig is strongly peaked around 
g = 0, as expected for a superfluid ground state, while it broadens as interactions are 
increased, as shown for u = 4. The shift ka is applied instantaneously, and the shape 
of Uq is unaltered at time t = O'^. This latter distribution is the initial condition for the 
subsequent analysis of the system dynamics. 

Figure [3] shows the current (j) calculated at t = 0+ for a few values of u. In the limit 
of vanishing interactions, the current (j) at t = O"*" is (j) = 2Jnsin [ka), corresponding 
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Figure 3. Boson current (j) in units where J = 1 after quasi-momenta < ka/ir < 0.5 
are imparted (time t ~ 0+). In the limit u ^ 0, (j) = 2Jsin(fca) (solid black line). 
Shown are results for interactions u = 1 and u — 4 with a filling factor of n = 1 
particles per site. The numerical parameters are x = 100, d = 6. 

to all particles occupying the same quasi-momentum state. In a MI state with integer 
n, this current will be zero in the limit of infinitely strong interactions, as the quasi- 
momentum distribution Ug becomes fiat. In figure |3] the case with m = is plotted as a 
solid black line. The figure shows that for finite values of the interaction strength u from 
figure [2], (j) decreases from the value 2 Jn sin (fca) as expected. However, we note that 
the current is comparatively large even for the initial condition u = 4, corresponding to 
a MI state at t = 0~. This is related to the fact that the quasi-momentum distribution 
Ug, whilst broader than in the superfiuid case, is still finite in width (see figure |2]), and 
far from the fiat distribution corresponding to the limit m — t- oo. As m is increased, the 
distribution continues to broaden gradually, but there is no strong qualitative change 
in the initial quasi-momentum distribution or the initial current values on entering the 
Mott Insulator regime. 

4-2. Time evolution of rig 

Due to the lack of Galilean invariance in a lattice, a state with a finite current does 
not correspond to an eigenstate of the system. We expect that in the presence of 
interactions any initial current (j) will eventually decay to a current-free state. The 
dynamics of the current decay corresponds to the rearrangement of the quasi-momentum 
distribution by two-atom scattering processes. The quasi-momentum of two particles is 
always conserved in these scattering processes, but because of Umklapp processes that 
connect the two edges of the Brillouin zone, the mean quasi-momentum is significantly 
changed, and will reach zero in a final steady state. 
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In figure S] we show results for tlie time evolution of Ug for a system with the initial 
conditions ka = OAn and u = n = 1. We observe for increasing time that the width 
of the momentum distribution increases and its mean value shifts towards zero. We 
have checked that the off-diagonal elements of the SPDM decay exponentially in the 
final state. These results are consistent with those reported in reference [23], where 
the change in the redistribution of quasi-momentum in the final state is linked to the 
appearance of a finite temperature in the system |33l HH SSI SS] . 

ka = OAtt 

0.05 
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q/lT 

Figure 4. Time evolution of the normalised ground state quasi-momentum 
distribution Uq in the SF regime for u = 1 after a momentum ka = OAn is imparted. 
The particles are redistributed into a final state with a broader distribution at smaller 
mean momentum on the timescale 10/J (see text). The numerical parameters are 
X = 100, rf = 6. 

We find the behaviour of and SPDM described above to be typical of all 
parameter values we have investigated, with the primary difference for different 
parameters being in the timescale on which these processes occur. This includes the 
Mott Insulator regime, where the main qualitative difference is in the timescale of the 
decay of the current, not in the form in which the quasi-momentum distribution is 
redistributed. The dependence of this timescale on the different system parameters will 
be investigated below. 

4-3. Stability diagram for the current 

We now investigate the decay of the current as a function of the initial mean quasi- 
momentum ka and the inter-particle interactions u for an initial ground state with 
n = 1. We focus on the time-dependent dynamics of (j) and also of the condensate 
fraction Cioo introduced in section 12.31 We extract a stability diagram for the crossover 
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between parameter regions where currents are stable over long time periods, and regimes 
where currents decay rapidly. 

In the following we focus on systems with interaction strengths 0.5 < u < 4.0 and 
O.OStt < ka < O.SOvr, and we compute the real time evolution of the boson current and 
the condensate fraction Cioo for times < t J < 10. This time interval is sufficiently 
short to allow currents to be probed in the experiment within typical decoherence times, 
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Figure 5. Time evolution of the boson current {j) for different initial parameters 
u and ka. In panel (a) results are compared for increasing values of the interaction 
u = 1, 2, 3, 4 for a constant small initial mean momentum ka = O.Itt. Panel (b) shows 
results for increasing initial momenta ka — O.Itt, 0.27r, O.Stt, 0.47r, O.Stt for u ^ 1. The 
numerical parameters are x = 100; d — Q. 
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e.g., due to incoherent light scattering [2], and also allows for accurate numerical results. 

Figure [5l^a) shows the time evolution of (j) for a constant initial momentum 
ka = O.lvr and increasing values of the on-site interaction u, and figure Et^b) shows 
the decay of (j) for a fixed on-site interaction of u = 1.0 and various values of ka. On 
the considered timescale we find that a constant superfiuid current (with no decay) only 
exists for small values of ka and u, for example ka = O.Itt with u = 1.0. For a slightly 
increased initial quasi-momentum or on-site interaction a small decay of about 10% on 
the timescale of the simulation becomes visible. Further increasing ka and/or u causes 
the decay rate to increase rapidly until for u > 3.0 at ka = O.Itt and ka > O.Svr at 
u = 0.1, the current has completely vanished at time tJ = 10. Figure [6] shows that 
the dependence of the time evolution of the condensate fraction Cioo on ka and u is 
qualitatively very similar to that of (j). 

An important result seen in both figure [5] and figure [6] is that we do not find a sharp 
transition point for u and ka separating parameter regions where the current is stable 
or unstable. This lack of a sharp transition is expected for one- dimensional systems, 
and is in contrast to the results found in higher dimensions [33] . 

Figures [7] and M summarise the dependence of the current decay on ka and u, and 
amount to stability diagrams for the crossover between parameter regimes where the 
current is stable and unstable. In particular, figure [7] is a shaded plot of the loss of the 
condensate fraction ACr{t), which we define as 

ACfi(r) = ICnit = r) - Cnit = 0)|/C^,(t = 0). 

We plot ACij(r) as a function of ka and -u at a specific time tJ = 10, as may be 
measured directly in an experiment. We chose R = 100, but found that ACr{t) is 
independent of R, based on a comparison of calculations with R = 50, 100, and 200 
sites. We see that stable currents exist on this timescale only for small values of u < 1.5 
and ka < O.lbTC, while no stable current is present for u > 2.5 and ka > 0.257r. In the 
intermediate region we observe a smooth crossover between those two regimes. In the 
classical limit of small on-site interactions, the stability/instability crossover tends to 
occur at large values of ka ~ 0.5'7r, which corresponds to the dynamical instability of 
reference [I9], while for fca ~ the instability sets in close to the value Uc ~ 3.37, which 
corresponds to the SF-MI transition at zero current [37]. As a reference, in figure [7] we 
plot the mean-field result of references [33l [3l] as a sohd black hne. The latter provides 
a reasonable indication of the position of the crossover region for small m < 0.5 only. 
Outside of this region, the decay appears typically to be much faster than is expected 
from the mean-field estimates. 

Note that ACji{t) discussed above can be directly observed via interference patterns 
in momentum distributions [1] , which can be measured via time-of-fiight measurements 
in experiments. We can estimate the corresponding experimental timescales by taking 
a typical lattice depth of 10 Er, where E^i denotes the recoil energy of the atoms 
Efi = p1/2m with the recoil momentum p^ = h/\i. For these lattice depths, the 
tunnelling amplitude J/h for Rb atoms is of the order of 100 Hz. Thus, the timescale 
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Figure 6. Time evolution of the condensate fraction Cioo evaluated from the SPDM 
over a range of 100 sites from the infinite homogeneous system. Qualitatively the 
behaviour is analogous to the boson current in figure [5] Panel (a) shows the decay 
for constant initial mean momentum ka and varying interactions u, and panel (b) for 
constant u and various ka. The numerical parameters are x = 100, d = 6. 



tJ = 10 corresponds to experimental timescales of the order of 100 ms, which is within 
typical coherence times. 

Figure [8] shows a time-independent stability diagram for the current, which we 
compute from decay rates Tc for the condensate fraction Cioo- As was shown in figure El 
the decay of Cioo as a function of tJ can be separated into three regions. Initially, 
the system shows an initial decay behaviour on short timescales tJ < 1, followed by 
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Figure 7. Shaded plot of the loss of the condensate fraction ACioo as a function of 
u and ka at a fixed time r J — 10. Values are shown on a u-ka grid with spacings 
of Au = 0.5 and Aka = 0.057r, smoothed by a spline interpolation. The mean-field 
Gutzwiller prediction from |33[I34[ for the transition is drawn as a solid black line. The 
figure visualises a stability crossover diagram (see text). The numerical parameters are 
X = 100, d = 6. 



an intermediate region where the decay is found to be approximately hnear in time. 
Finally a saturation of the decay process takes place while the system approaches the 
zero-current steady state. We therefore extract the decay rate Tc by determining the 
slope of the decay in the approximately linear intermediate region. Note that in the short 
time region where we fit, a linear decay behaviour is also equivalent to an exponential 
decay Cioo(^) oc exp(-rct). 

Figure [HI which is a shaded plot of F^, for F^/J < 0.1, shows results which are 
qualitatively similar to those of figure [71 That is, stable currents are found to exist only 
for small values of ka and u, and the stability /instability crossover in general occurs at 
values of ka and u significantly smaller than predicted by mean-field theory (see solid 
black line in the figure). 

In the next section we compute decay rates for various lattice fillings n > 1 
and investigate the scaling with the interaction strength and filling factor. We find 
suprisingly that the scaling is very similar to those computed for a weakly interacting 
system with fi 3> 1 using beyond-mean-field (phase-slip) calculations. 

5. Scaling of the current decay rates 

In this section, we study the dependence of the decay rate Tq defined in section 14.31 
and of an analogous time-independent decay rate Tj for the current (defined below) on 
the lattice filling n, the interaction strength u, and the quasi-momentum ka. We find 
scaling laws exhibiting exponential dependence on these quantities. Interestingly, these 
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Figure 8. Shaded plot of the condensate fraction decay rates Tc/J as a function of 
u and ka. The decay rates are extracted from the time evolution of the condensate 
fraction Cioo by a linear fitting technique (see text for details). Values are shown 
on a u-ka grid with spacings of Au — 0.5 and Aka = O.OStt, smoothed by a spline 
interpolation. All decay rates Tc/J > 0.1 are marked in white colour. The Gutzwiller 
mean-field prediction is drawn as a solid black line. The figure is in qualitative 
agreement with figure [7] and visualises a time- independent stability crossover diagram 
(see text). The numerical parameters are x ~ WO, d — 6. 



scaling laws are similar to those computed for current decay in a weakly interacting 
system with large filling factor n ^ 1, which were computed via (beyond- mean-field) 
instanton calculations in reference [M]- This is especially true in the case n > 3, 
or far from the SF-MI transition with n = 1. Naturally, significant deviations occur 
between the exact numerical results close to the SF-MI transition and predictions 
for a weakly interacting system. Note that here we do not attempt to verify these 
instanton predictions within their regime of validity, which has been done in a recent 
study by Danshita and Polkovnikov [17], using particle numbers n ~ 1000. Instead, 
we are interested in exploring the regimes close to n = 1, which correspond to recent 
experiments with atoms confined in 3D optical lattices. 

In reference [M], a decay rate F for the current in a system with large filling factor 
was calculated using instanton (phase-slip) techniques. This could be shown to be 
F oc e~'^ , with 5* the semi-classical action 



Equation (jH]) is strictly valid in ID in the weakly interacting (Gross-Pitaevskii) limit, 
u/fi <^ 1, for un ^ 1, i.e. in the limit of large filling factors, and close to the classical 
instability at ka ^ tt/2. This is very different from the parameter regimes that are 
typical when the ID system is formed via a 3D optical lattice, as we study here in 
the present article. In principle, our filling factors n < A, and acceleration values 
0.2 < ka/n < 0.35, do not correpond to the region of validity of equation (jG]). However, 




(6) 
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as we will see below, we nonetheless obtain very similar scalings. 
6 
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Figure 9. The negative logarithm of the decay rate of Cioo in panel (a), and of the 
boson current in panel (b) as a function of \ffiju for a fixed initial mean momentum 
ka = 7r/4. The results are estimated from linear decay rate fits. Both panels visualise 
a similar linear dependence (see text). The numerical parameters are x = 50, d = 10. 



Figures [9](a) and[9](b) show our results for the decay rates Tc and Tj, respectively, 
where Tj has been extracted from the time evolution of the boson current (j) (see 
figure [5]) in the same way as for Tq as described in section 14.31 In panels (a) and (b) 
we plot — ln(rc/J) and — ln(rj/J) respectively as a function of ^/nju^ for 1 < n < 4 
and a fixed ka = n/A. Both of these plots show an approximately linear dependence on 
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y^n/u, for all n and \ffiju > 0.6. In addition, for these values of \/n/u both panels 
show that by increasing n the slope of the lines tend to converge to a fixed value. This 
indicates that not only is the decay rate proportional to exp(— A-^I/m), but that the 
constant A appears to scale as y/ft for n > 2. This is the same scaling as predicted for 
the case 1 via instanton methods. Averaging over the slopes of the results for the 
large filling factors n = 3, 3.5 and 4, in order to obtain the prefactor in A, we find values 
of A = (2.1 ± 0.5)a/^ and A = (2.7 ± 0.6)a/^ for panels (a) and (b) respectively. We 
note that the equality of the slopes obtained from the current decay and condensate 
fraction decay within the fitting error indicates the equivalence of these quantities as 
good dynamical observables, for the characterisation of the current decay both in the 
numerical simulations and in experiments, for large enough n. For small n ^ 1 we 
observe slightly different behaviour (see below for more details). 

The results above for large n surprisingly show good agreement with scalings from 
the instanton calculations for a very different parameter regime, although we note that 
the prefactor in A is different, as the scaling law in equation (E]), would imply A = 3.9\/^- 
The most obvious deviation between the iTEBD results and equation ([6]) takes place 
for the Tq results of panel (a) in the vicinity of the SF-MI transition with n = 1, 
which occurs at \ffifu ~ 0.54. For \ffifu < 0.6 the rate no longer increases as 
\/n/u decreases. We indicate this parameter region with an arrow in figure [9](a). This 
behaviour is not captured by the decay of the current F^ in panel (b). This contrast 
between the the decay of the condensate fraction and the decay of the current could 
refiect a reduced role of the condensate fraction in the overall system dynamics when 
we have an initial MI state. Note that for \ffifu < 0.6 the decay of (j) occurs rapidly 
on the timescale of our calculation time steps, which makes computing its decay rate 
less accurate in this regime. 

6. Summary and Outlook 

In summary, the iTEBD method allows us to make quantitative predictions for the 
time- dependence of currents for bosons moving in ID in an optical lattice. In contrast 
to the behaviour in higher dimensions, we observe a broad crossover between regions of 
stable and unstable currents, with the typical decay rates for the currents being more 
rapid than mean-field predictions. We also find surprising agreement with the scaling of 
decay rates with interaction strength and filling fraction calculated for the case of large 
fillings and weak interactions via phase-slip methods. These results should be accessible 
in current experiments with filling factors near unit density. As the density increases, 
we find that the scaling of the decay rate of the current almost agrees with predictions 
from instanton calculations. 

These results are strongly related to the decay of currents generated in a displaced 
harmonic trap, which has also been the subject of significant experimental [161 113 EO] 
and theoretical [191 ISH [22] investigation. An understanding of these currents is also 
an important starting point for the investigation of more general transport properties 
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of bosons in ID. This could include behaviour of these currents in the presence of an 
impurity [3 [8], or spatially varying interactions |9l [TO] . 
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Appendix A. Translationally invariant block size and the iTEBD Algorithm 

The original iTEBD algorithm as introduced in [35] uses translationally invariant matrix 
product states to represent the quantum state of an infinite homogenous system. In order 
to conveniently compute time-evolution, the translationally invariant state is usually 
represented by a repeated block of two sites, as described in reference [35] . 
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Figure Al. The time evolution of Cioo after an initial quasi- momentum acceleration 
ka ~ 7r/4 for two different on-site interactions u = 1 and w = 3. No differences are 
visible when doubling or tripleing the block size m within the iTEBD algoritm. Other 
numerical parameters are x — 100, d = 6. 

Note that despite the block size of two, excitations with any period can be 
represented in this form, because matrix product states effectively store phase 
relationships between neighbouring sites rather than phases corresponding to a single 
site. Indeed, the use of a block size of two is a matter of convenience in applying the 
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algorithm, and if variational methods are applied [TT], then a block size of one can be 
sufficient to represent a translationally invariant state. In the present work, we find 
it convenient to extend the block size beyond two to m sites, where for a particular 
simulation in which the momentum translation of the initial state is ka, we choose 
m = 27r/ ka. This is convenient in the application of the operator K{ka), but has no 
fundamental effects on the computation of time evolution of the state. In our case, 
time evolution is simulated by decomposing the time evolution operator into m two- 
site unitary operations, instead of two as in the original iTEBD algorithm. We expect 
that this, in case of a large enough bond dimension x still results in an exact state 
representation for the infinite system during the whole time evolution. 

To demonstrate the independence of the computation of time evolution on our 
choice of m, we plot example values for the time evolution of the condensate fraction Cioo, 
computed with increasing m in figure lATl We use an initial acceleration ka = 7r/4, and 
compare results for block sizes m = 8, 16, 24, for two different values of the interaction 
strength u = 1,3. We see that the results are independent of the block size, and after 
testing for convergence in the parametes x and d should represent exact time-evolved 
values. 
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